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Abstract 



Estimation of division and death rates of lymphocytes in different conditions is vital for quan- 
titative understanding of the immune system. Deuterium, in the form of deuterated glucose 
or heavy water, can be used to measure rates of proliferation and death of lymphocytes in 
vivo. Inferring these rates from labeling and dclabcling curves has been subject to considerable 
debate with different groups suggesting different mathematical models for that purpose. We 
show that the three models that are most commonly used are in fact mathematically identical 
and differ only in their interpretation of the estimated parameters. By extending these previous 
models, we here propose a more mechanistic approach for the analysis of data from deuterium 
labeling experiments. We construct a model of "kinetic heterogeneity" in which the total cell 
population consists of many sub-populations with different rates of cell turnover. In this model, 
for a given distribution of the rates of turnover, the predicted fraction of labeled DNA accumu- 
lated and lost can be calculated. Our model reproduces several previously made experimental 
observations, such as a negative correlation between the length of the labeling period and the 
rate at which labeled DNA is lost after label cessation. We demonstrate the reliability of the 
new explicit kinetic heterogeneity model by applying it to artificially generated datascts, and 
illustrate its usefulness by fitting experimental data. In contrast to previous models, the ex- 
plicit kinetic heterogeneity model 1) provides a mechanistic way of interpreting labeling data; 
2) allows for a non-exponential loss of labeled cells during delabeling, and 3) can be used to 
describe data with variable labeling length. 



Classification: Immunology, Mathematical Biology. 

Keywords: D-glucose, parameter estimation, cell turnover, lymphocyte dynamics, heavy wa- 
ter 

Short title: Measuring cell turnover 

Abbreviations: D-glucose = ^i72-glucose, heavy water = ^H20, KH = kinetic heterogeneity, 
AM = asymptote model, RSS = residual sum of squares. 



1 Introduction 



There is little consensus about the expected life spans of lymphocyte populations in health 
and disease. Labeling the DNA of dividing cells with deuterium has proved to be one of the 
most reliable and feasible ways to study the population dynamics of lymphocytes in healthy 
human volunteers and in patients (jH, 12 S)- Deuterium, in the form of deuterated glucose or 
heavy water, is used to measure the rate at which cells are dividing in vivo, without the need to 
interfere with these cellular kinetics. Deuterium is incorporated into newly synthesized DNA 
via the de novo pathway (0), and enrichment of deuterium (over hydrogen) in the DNA of cells is 
therefore related to cell division. During label administration, the fraction of deuterium-labeled 
nucleotides increases over time, and after label withdrawal, the fraction generally declines over 
time (0, 0). Labeling DNA with deuterium in humans has a number of clear advantages 
over other labeling techniques such as with BrdU, including the absence of toxicity, the fact 
that the rate of incorporation of deuterium into the DNA is independent of the amount of 
nucleotides present, and a simpler mathematical interpretation of the data (0, [s], 0). Several 
mathematical models have been proposed for estimation of cellular turnover rates from labeling 
data (reviewed in (0)), including a simple precursor-product relationship (El), a source model 
(0), a two-compartment model (jl, [o]), and a kinetic heterogeneity model (jiol ). 



In their pioneering study on deuterium labeling, Mohri et al. (|2|) found that the estimated 
rate of cell proliferation was typically smaller than the rate of cell death. Because the cell 
population under investigation was in steady state, the extra death must be compensated by a 
source of cells, for example from the thymus. This interpretation was challenged by the work of 



Asquith et al. (llOl ). which pointed out that estimated proliferation and death rates do not have 
to be equal if the population is kinetically heterogeneous (i.e., different cells in the population 
divide and die at different rates). Because the labeled population preferentially contains cells 
that proliferate (and die) relatively rapidly, the estimated rate of cell death is in fact expected 
to be higher than the average proliferation rate (jlol ). 

Here we extend these studies and propose a more mechanistic approach to estimate the 
rates of lymphocyte proliferation and death from deuterium labeling experiments. First, we 
show that the three most commonly used mathematical models are in fact kinetically identical 
(i.e., will lead to identical estimates of the average rate of cell turnover), and only differ in their 
biological interpretation of the model parameters. Second, we formulate a novel mathematical 
model which explicitly takes into account kinetic heterogeneity of lymphocyte populations, and 
show how lymphocyte turnover rates can be calculated using this model. Several previously 
made experimental observations arise naturally from the new model. For example, we find 
that the rate of label loss during delabeling generally exceeds the rate of label accumulation 
during the labeling phase. Our model also explains the dependence of the rate at which labeled 
DNA is lost after label withdrawal on the duration of the labeling period (jToh . As a proof of 
principle, we demonstrate that the newly developed model can fit artificially generated data. 
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and correctly returns their underlying kinetic parameters. We also illustrate the usefulness 
of the new model by fitting it to several experimental datasets. The novel explicit kinetic 
heterogeneity model may offer alternative interpretations of how infections or treatments affect 
the turnover of human lymphocytes in vivo. 



2 Results 



2.1 Previous models 



Although different models have been proposed for interpretation of deuterium labeling data 
(0, [lol ) and are being debated in the literature, they are in fact mathematically identical. 
Following De Boer et al. (jlll ). consider a cell population consisting of a fraction a of cells with 
average turnover rate d (i.e., an expected life span of l/d days), and a fraction 1 — a of cells 
that do not turnover at all on the time scale of the experiment. During the labeling phase, 
consider the fraction of unlabeled DNA [/„ in the sub-population with death rate d. Because 
DNA is only lost by cell death, [/„ changes according to: 



df/„/dt = -dUa. 



During the delabeling phase the fraction of labeled DNA in that same population (Lq,) is 
described by: 

dLa/dt = —dLa, 



because labeled DNA can only be lost by cell death. Since Ua + = 1, the fraction of labeled 
DNA in the whole population L{t) = aL^it) is described by: 



r «(l-e-'^*), ift<T, 
^^^> ~ \ L(T)e-'^(*-^), otherwise, ^ ' 

where T is the duration of the labeling period. Given that only a fraction a of all cells in 
the population are turning over (or dying) at rate d, the average turnover rate of the whole 
population is a x c? (Tl'). Importantly, this approach does not require us to describe how new 
cells are formed, i.e., they could be generated by the thymus and/or by proliferation. Extending 
this model by assuming n sub-populations with different rates of cell proliferation pi and death 
di, and possibly generation of new cells from a source Si (Figure [1]), the fraction of labeled 
nucleotides in the whole population at time t is given by: 
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^ ' ~ I Eti o^i (1 - e-*"") e-^-(*-^), otherwise, 

where is the fraction of cells in population i with death rate di, and a = J2i=i < 1 is the 
asymptote that would be approached if label would be administered indefinitely. 

The "source" model that was previously proposed by Mohri et al. (i2i) considered one homo- 
geneous cell population, but allowed for a source of unlabeled cells during the labeling phase, 
i.e., dU/dt = su — dU, which also gives rise to an asymptote a = 1 — su/d, defining the fraction 
of cells that can maximally become labeled (here su = su / X where X is the total number of 
cells in the population at equilibrium and su is the number of cells with unlabeled DNA coming 
from the source per day during the labeling phase (ji])). Mathematically, the source model is 
therefore identical to eqn. ([T]). Similarly, in the kinetic heterogeneity model devised by Asquith 
et al. (Iiol ). dL/dt = p{U + L) — dL = p — dL for the labeling phase and dL/dt = —dL for 
the delabeling phase. Assuming p < d one again obtains eqn. ([T]) with a = p/d. Therefore, 
all these models are mathematically identical and only differ in the biological interpretation of 
their model parameters (see also (]l2|)). We propose to call all these models the "asymptote 
model". Importantly, in all models the product a x d can be interpreted as the average rate 



of cell turnover of the population as a whole (jlll ). and therefore, all three models, when fitted 



to data, will deliver identical estimates of the average turnover rate, which is the parameter of 
interest . 



2.2 Kinetic heterogeneity model with continuously distributed turnover 
rates 



Because of its simplicity, the model given in eqn. ([T]) has two limitations. First, the as ymp tote 
level is a phenomenological parameter that depends on the length of the labeling period (llOl ). As 
a consequence, datasets with different labeling periods will likely give rise to different estimated 
asymptotes and different estimated average rates of cell turnover. Therefore strictly speaking, 
this model cannot be used to explain multiple datasets coming from the same experimental 
setup varying only in the length of the labeling period; the differences in the rate at which 
labeled DNA is lost would force either the asymptote or the estimated average turnover rate to 
be different for the different labeling periods (Den Braber et al. in preparation). Second, the 
model assumes that the increase in labeled DNA during the uplabeling phase, and the loss of 
labeled DNA during the delabeling phase can be described by single exponentials. This may 
be incorrect if cell populations with different turnover rates are labeled and subsequently lost. 

Under very general assumptions, we have formulated an alternative model that does not 
make these a priori assumptions. In our new model, a cell population consists of n sub- 
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populations each with different kinetic properties (see Figure [U and eqn. ([2])). If the number 
of sub-populations is large (n — *• oo), the sum in eqn. ([2]) can be replaced by an integral. The 
fraction of labeled nucleotides in the population then becomes (see Supplemental Information 
for derivation) 



L{t) = / L{t,d)dd (3) 
Jo 

where L{t,d) is given by eqn. ([1]) where a = f{d) is the frequency distribution of turnover 
rates, and f{d)dd is the probability that a randomly chosen cell in the population belongs to 
a sub-population with a turnover rate in the range {d,d + dd). If the turnover rates in the 
population, f{d), follow a gamma distribution, the change in the fraction of labeled DNA with 
time is given by: 



dt 



1 + f) , if t < T, 

1 + ) - ( 1 + f ) ' otherwise. 



where d is the average rate of cell turnover in the population, k is the shape parameter of 
the gamma distribution, and T is the duration of the labeling period. For k = 1, the gamma 
distribution becomes an exponential distribution, and the rate at which the fraction of labeled 
DNA changes is simply: 



dt 



if t < T, 



L{t) = { 1+'^*' ™ - (5) 

(i+,,)(i+rf(,-T)) ^ Otherwise. 

This is an interesting model in which a single parameter J predicts both the rate of uplabeling 
and downlabeling, and in which there is no asymptote below 100% for the level of labeled 
DNA, i.e., under continuous label administration all cells in the population will become labeled 
(Figure [2]). Moreover, this model predicts that the initial rate d* at which labeled DNA is lost 
after label cessation depends on the duration of the labeling period, d* ^ d{l + [1 + dT)~^) (see 
Supplementary Information for derivation). According to this model, short labeling experiments 
{dT <^ 1, d* ^ 2d) will lead to 2- fold faster initial rates of decline in the fraction of labeled 
nucleotides than longer labeling experiments {dT ^ 1, rf* ^ rf). 

Solutions (jlj) and ([5]) predict that the initial rate of increase in the fraction of labeled DNA is 
the average rate of cell turnover d (see also Supplementary Information). However, the increase 
in the fraction of labeled DNA does not appear to be exponential, as was implicitly assumed 
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in the asymptote models discussed above. Similarly, during the delabeling period, the model 
predicts a non-exponential decline in the fraction of labeled DNA (Figure [2] and [3]). In general, 
the initial rate of label loss during delabeling d* is given by: 



d* ^d + 



var((i) 



(6) 



d 



where \ai{d) is the variance of the distribution of turnover rates in the population. In case 
when turnover rates follow a gamma distribution, the initial rate of loss of the label after short 
labeling periods depends on the shape parameter k of the distribution, d* d{l + 1/k), while 
it does not after long labeling periods [d* ~ d). The rate of loss of labeled DNA slows down as 
less DNA remains labeled, which is most clearly seen when proliferation rates are distributed 
according to a very skewed gamma distribution (A; < 1, Figure [3]). This is a natural property 
of the explicit kinetic heterogeneity model as loss of labeled DNA is reflecting the distribution 
of the turnover rates of the different sub-populations, with labeled DNA from the most rapidly 
turning over sub-populations being lost first (early fast decline) and labeled DNA from the 
other, more slowly turning over, populations being lost later (late slow decline). 

To study the effect of the shape of the turnover rate distribution on the predicted labeling 
curve, we plotted the changes in the fraction of labeled DNA as predicted by the model (Figure 
[3JA.&B) with different gamma-distributed turnover rates (Figure [3]C). When the gamma distri- 
bution is highly skewed (i.e., A; < 1), the majority of cell sub-populations have very low rates 
of cell turnover, and the average rate of cell turnover is dominated by a few sub-populations 
that turn over unrealistically fast. This is best illustrated by calculating the cumulative con- 
tribution of a sub-population with a particular rate of turnover to the average turnover rate of 
the population d (Figure [3p): 



For large values of the shape parameter k (e.g., = 5), the sub-populations with turnover 
rates that are close to the average turnover rate are the main contributors to the average rate 
of cell turnover (Figure [3p). When the gamma distribution is extremely skewed [k = 0.01), 
the rate of turnover of the sub-populations that contribute significantly to the average turnover 
rate is as high as 10^ — 10^ per day, which is biologically unrealistic. Therefore, the gamma 
distribution should be rejected whenever one estimates a high average turnover rate d with 
a highly skewed gamma distribution (i.e., a low value of the shape parameter k). As a rule 
of thumb, k should be larger than 0.1 (Figure [3p); otherwise a relatively large fraction of 
sub-populations has unrealistically high turnover rates. 




(7) 
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It is possible, however, that not all cells in the population are turning over. The models 
above can easily be extended to incorporate this possibility by allowing for the same asymptote 
as in eqn. ([1]). An example would be a labeling experiment in which slowly turning over naive 
T lymphocytes and more rapidly turning over memory lymphocytes are not separated (0). If 
only a fraction a of cells have turnover rates that are distributed exponentially, and the other 
cells undergo negligible turnover on the time scale of the experiment, the change of the fraction 
of labeled nucleotides with time is given by: 



adgt \f f < T 

L{t) = { 1+"'^*',-^ , ' (8) 

\ J > ""^-^ otherwise, ^ ^ 



(l + dc,t){l + do,{t-T))' 



where da is the average of the exponentially distributed turnover rates, and the average rate of 
cell turnover in the whole population is d = ad^- 

It should be noted that the results of this section are applicable both to proliferating and 
non-proliferating lymphocytes, given the general structure of the cell population in the model 
(see Figured]). As a downside of this, the model does not allow to estimate which fraction of 
labeling of lymphocytes is due to proliferation of precursors (e.g., thymocytes for naive T cells) 
or due to peripheral proliferation of the lymphocyte population itself. Additional experiments, 
such as thymectomy in case of studies of naive T-cell turnover, may allow to estimate the 



separate contribution of peripheral T-cell proliferation (jl3l . Den Braber et al. submitted) 



2.3 Fitting artificial data to validate the model 

Having analytical expressions for several kinetic heterogeneity models, we analyzed how well 
these models can recover the (known) average turnover parameter from simulated (artificial) 
datasets. Three models were used to generate artificial datasets: 1) the kinetic heterogeneity 
model with gamma-distributed rates of turnover (eqn. (j3])), referred to as the "Gamma model", 
2) the kinetic heterogeneity model in which a fraction a < 1 of cells have exponentially- 
distributed rates of turnover (eqn. ([8])), referred to as the "Exponential model", and 3) a 
"Two population model" (eqn. ([2]) with n = 2, turnover rates di and d2, average turnover rate 
d = adi + (1 — (y)d2, and di > d). These datasets were subsequently fitted by the same three 
models as well as by the conventional Asymptote model (eqn. ([T])). 

Not surprisingly, the models delivered correct estimates for the average turnover parameter 
if a dataset was fitted with the model that was used to generate the data (Tableland Figure [5l 
with the results obtained by fitting the Two population model not shown). All models described 
the data sets generated by the other models reasonably well (Figure H]), although some features 
in the data could not be reproduced. For example, the Asymptote model failed to describe the 
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decreasing rate at which labeled DNA is lost over time, which is observed in the data generated 
by the Gamma and the Exponential models (see last data points in Figure H]). Some model 
fits delivered incorrect estimates for the average turnover rate if the data were generated using 
another model. For example, the Gamma model overestimated the average turnover rate when 
the data were generated using the Exponential model (up to 2-fold), and underestimated d 
for data generated using the Two populations model (over 2- fold). This is most likely due 
to the strong constraint of the model that both uplabeling and delabeling curves have to be 
described with one mechanism, i.e., gamma-distributed turnover rates. On the other hand, the 
Asymptote model always underestimated the true average turnover rate (up to 2-fold for data 
generated by the Two-populations model; Tabled]). It did perform somewhat better than the 
Gamma model as judged by the mean square distances, because the rate of uplabeling and 
downlabeling are relatively independent in the Asymptote model. 

Given that natural lymphocyte populations are likely to contain resting sub-populations, 
some extent of saturation in the fraction of deuterium-labeled nucleotides is expected in almost 
any experimental dataset. In our artificial data, such an asymptote was imposed when using 
the Exponential model by letting only 50% of all cells to turn over (Figure [Hand Tabled]). It is 
therefore not surprising that the Gamma model, which does not have an explicit asymptote in 
the uplabeling phase (see eqn. (jl])), did not correctly estimate the average turnover rate for the 
data generated by the Exponential model (Figure [5]). Extending the Gamma model to allow 
for an explicit asymptote during the labeling phase [a) indeed improved the estimate of the 
average turnover rate d = ada = 0.13 day~^ (with 95% CIs = 0.095 — 0.22 day~^ which includes 
the true average d = 0.1 day~^), even though the estimated fraction of turning over cells a was 
not significantly different from 1 (i.e., an F-test would not reject a model with a = 1; results 
not shown). This exercise illustrates that when fitting experimental data one should check 
whether allowing for an explicit asymptote in the uplabeling phase leads to different estimates 
of the average turnover rate. Interestingly, all models underestimated the average rate of cell 
turnover when the data were generated using the Two populations model. This is because the 
models did not reproduce the relatively rapid accumulation of the labeled DNA in the first 
days (Figure H]). Fitting the Two populations model to these data led to better estimates of the 
average turnover rate {d = 0.084 per day with 95% CIs = (0.062 — 0.35) for 7 days of labeling, 
and d = 0.26 per day with 95% CIs = (0.061 — 0.44) for 15 days of labeling, where the constant 
d = 0.1 day~^ is contained within both ranges, results not shown). 

Although stable isotope labeling seems to be the best tool at hand to estimate rates of 
lymphocyte turnover, a recent review (0) pointed out that estimated lymphocyte turnover 
rates differ consistently, depending on the labeling method used (heavy water or deuterated 
glucose), and the length of the labeling period. A priori., according to the Asymptote model 
that is generally used, the estimated average turnover rate should not depend on the length of 
the labeling period. Using our explicit kinetic heterogeneity models, we analyzed the influence 
of the length of the labeling period on the estimated average turnover rate. For all models, 
we found that the duration of labeling had little influence on the estimated average turnover 
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rate (for the chosen labehng periods of 7 and 15 days, Table [1]). This suggests that longer 
labeling periods will not necessarily result in lower estimates of the average cell turnover rate 
than shorter labeling periods. 

An overall conclusion of this analysis is that without a good understanding of the underlying 
model of cell proliferation (i.e., the distribution of turnover rates in the population), one may 
obtain incorrect estimates of cellular turnover rates, even if the quality of the fit to the data 
is acceptably good. Therefore, when analyzing experimental data, one should aim at using 
several alternative models for fitting, and investigate whether estimates of kinetically important 
parameters, such as the average rate of cell turnover, are independent of the model used. 



2.4 Fitting experimental data 



We next sought to determine how well the new kinetic heterogeneity models fit experimental 
data. Using deuterated glucose, Mohri et al. ^ obtained labeling data of T lymphocytes from 
uninfected healthy human volunteers and from chronically HIV-infected patients. Previously, 
these data were fitted using an extended 4-parameter source model, to estimate the rates of 
cell division and death of T lymphocytes in healthy humans, and to obtain insights into how 
these rates change upon HIV- infection (2). Lymphocytes were sorted into CD4+ and CD8"*" T 
cells, without distinguishing between their naive and memory subpopulations. Since naive T 
cells have a much slower rate of turnover than memory T cells (jsl), it is natural to assume an 
asymptote in the fraction of labeled nucleotides of unsorted CD4+ and CDS"'" T cells. 



We have refitted the labeling data from the four healthy human volunteers studied by Mohri 
et al. (2), again using the three models for cell proliferation: the Asymptote model (eqn. ([1])), 
the Exponential model (with a fraction a of cells with exponentially distributed turnover rates, 
eqn. ([8])) and the Gamma model (with gamma distributed turnover rates, eqn. (jlj)). The 
data were fitted simultaneously for all four healthy volunteers while searching for the minimal 
number of parameters that describe the data with reasonable quality (using a partial F-test 
for nested models (Il4l)). Overall, the models described the data reasonably well (Figure [6] and 
Table [2] and [3]). For CD4+ T cells, the average turnover rate and the delay at which labeled cells 
appeared in the blood did not differ significantly between patients (Table [2]). For all patients, 
the average rate of turnover was about 0.46% per day with a corresponding estimated half-life 
of In 2/ J ^ 151 days. There was an average delay of one day before labeled cells appeared 
in the blood. The average turnover rate of 004"*" T cells from control cl was always higher 
than that of the other individuals, irrespective of the model used (Figure [7]/\.), which may be a 
sign of an immune response to an infection in cl (see also below). Both the Asymptote model 
(a ~ 15%) and the Exponential model (a ~ 25%) predicted an asymptote in labehng that is 
smaller than the fraction of memory phenotype 004"*" T cells in humans of that age (}3|). The 
Gamma model could describe these data even better than the other two models with no need 
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for an asymptote. The observation that the estimate of the asymptote can differ dramatically 
between different models reconfirms our statement that this parameter is of little use for data 
interpretation (jlll ). 



For CD8"'" T cells, the parameters differed significantly between different healthy volunteers, 
with the exception of the asymptote level a in the Exponential model which could be fixed 
between individuals. The estimates of the average turnover rates of CDS"*" T cells in healthy 
volunteers c2-c4 did not strongly depend on the model that was used to fit the data. However, 
the estimated turnover rate of CD8"'" T cells in individual cl, which was much higher than 
the estimated turnover rate in the other healthy volunteers, depended strongly on the model 
used and was estimated to be the highest when using the Gamma model. The latter model 
fitted the labeling data from all four individuals very well and reproduced the non-exponential 
change in the fraction of labeled DNA in the downlabeling phase (Figure EF). In all four healthy 
individuals, CDS"*" T cells turned over at a slower rate than CD4"'' T cells; the average turnover 
rate of CD8^ T cells was d = 0.29% per day with a corresponding half-life of ln2/(i ^ 239 days. 
The fits of the Asymptote model and the Exponential model predicted an asymptote in labeling 
of a ~ 0.15 (Table[2]). Even though the Gamma model lacks an explicit asymptote lower than 1, 
it fitted these data with equally good quahty as the models with explicit asymptotes. Allowing 
for an explicit asymptote in the Gamma model did not improve the quality of the fit (CD4"'" 
T cells: p > 0.38, CD8+ T cells: p > 0.99, F-test), and the estimated average lymphocyte 
turnover rates were not affected by the addition of an explicit asymptote (results not shown). 



It is important to investigate whether the good description of the data of the model with 
gamma distributed turnover rates is achieved with biologically reasonable parameter values. In 
all data we estimated the shape parameter of the gamma distribution to be small, i.e. k < 1, 
but k was estimated to be larger than 0.1 in seven of the eight fits. Low values of the shape 
parameter k imply that in the population most cells turn over at very slow rates while a few 
populations turn over very rapidly. To investigate whether such a distribution is biologically 
reasonable, we calculated the fraction of cells in the population with a turnover rate higher than 
dmn^r = 1 per day which is the maximal rate of CD8+ T-cell proliferation in rhesus macaques 
(llSl). This fraction is given by f{d)dd for the estimated parameters of the distribution 
(see Table [2] and [3]) . For most fits, the fraction of cells with turnover rates higher than 1 per 
day is < 10~^^, and given the estimated total number of lymphocytes in humans of ~ 10^^ (jl6l ). 
that would yield only a few cells with unrealistically high rates of turnover. However, for the 
CD8+ T cells of healthy volunteer cl we found that ~ 3 X 10"^ X 10^2 = 3 X 10^ cells turn over 
at rates higher than 1 per day, which is unrealistically high. 



To investigate this further we reanalyzed the CD8+ T-cell labeling data of individual cl 
using several extended models. In the first model, a fraction a of cells in the population have 
gamma-distributed turnover rates while the other fraction (1 — a) of cells turn over at the 
highest possible rate dmax = 1 day~^. This situation may correspond to a scenario where a 
small fraction of CD8+ T cells is responding to an infection. However, this model failed to 
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describe the data with biologically reasonable parameter values (a ~ 1 and k = 0.03). 



In the second extended model, the gamma distribution of turnover rates was truncated at a 
maximal value dmax = 1 (see Supplementary Information for analytical results). The fit of this 
model to the labeling data for individual cl was of similar quality as the fit in which the gamma 
distribution was not truncated, and it delivered similar estimates for the average turnover rate 
and the shape parameter {d = 0.66% per day and k = 0.030, results not shown). We estimate 
that in healthy volunteer cl about 0.1% of all CD8^ T cells are turning over rapidly at rates 
between 0.5 — 1.0 per day, which is reasonable. For example, in mice responding to lymphocytic 
choriomeningitis virus (LCMV) infection, at the peak of the immune response more than 50% 



of all CD8^ T cells in the spleen are specific for the virus (llTl. llSl). 



Finally, in the third model, we assumed that the CD8^ T-cell population in patient cl 
consists of naive, memory and effector T-cell subpopulations with 3 different rates of turnover 
(see eqn. ([2]) with n = 3). Assuming that the rate of turnover of naive T cells is and that 
letting for effector cells dmax = 1 per day, we could obtain excellent fits of the labeling data with 
an estimated average turnover rate d = 1.4% per day (95% CIs = (1.2 — 1.6)% per day) which 
is much higher than estimates obtained by other models (Figure [7]). Using model selection 
methods such as the Akaike Information Criterion, we found equal support for the latter model 



and the model in which the turnover rates follow a truncated gamma distribution (Il9l . results 
not shown). We can conclude, therefore, that the average turnover rate of CDS"'" T cells in 
patient cl is at least 0.62% per day (Gamma model) and could be as high as 1.4% per day 
(Three population model). In summary, it seems that the average turnover rate of both CD4"'" 
and CDS"'" T cells was increased in individual cl as compared to other individuals, and this 
could be explained by a normal immune response in this otherwise healthy volunteer. 



3 Discussion 



In this paper we have analyzed the models that are commonly used in the literature to estimate 
the rates of cell turnover from deuterium labeling data. We have shown that the three most 
commonly used models are mathematically identical and therefore provide identical fits to the 
data, and only differ in the biological interpretation of the estimated parameters (jl2l ). The 
simplest summary of labeling data is provided by a model that has two parameters: d as the 
rate of cell death in the population, and a as the fraction of cells that undergo turnover, which 
determines the asymptote of the uplabeling phase (see eqn. ([1])). In this model, a x d gives the 
estimated average rate of cell turnover. We have extended this model by allowing for multiple 
sub-populations i of size with different turnover rates di (see eqn. ([2])). This extended 
model can be used to investigate potential heterogeneity of cell populations, by fitting labeling 
data with a model that has one, two, or more sub-populations with different turnover rates. 
Using standard techniques of model selection (e.g., the partial F-test or the Akaike Information 
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Criterion), one can investigate which of those models describes the labehng data best, given the 
number of model parameters (1l4l . Il9l). or one can study whether the estimated average turnover 
rate is converging to an invariant value by increasing the number of compartments (work in 
progress) . 



For the case where the number of sub-populations is large, we derived a model with con- 
tinuous kinetic heterogeneity. For several continuous distributions such as the exponential and 
the gamma distribution, the model predicts that the initial rate of loss of labeled DNA after 
label withdrawal is determined by the duration of the labeling period as has been observed 
experimentally (iq|)- Moreover, in the model the average turnover rate, which determines the 
initial rate of label accumulation in the population, turned out to be independent of the length 
of the labeling period. However, it should be noted that the average rate of cell turnover that 
is estimated from experimental data using, for example, the Asymptote model, may in fact 
depend on the duration of labeling (0, Den Braber et al. (in prepartion)). Potential reasons 
for this discrepancy will be investigated in more detail elsewhere. 

Previous models had certain artifacts: the asymptote labeling level was dependent on the 
length of the labeling period, and the accrual of labeled DNA during the uplabeling phase 
and the loss of labeled DNA during the downlabeling period were always described by single 
exponential functions. It is, therefore, unclear whether such limited models provide a good 
description of truly kinetically heterogeneous populations. Although at first glance it may seem 
impossible to fit labeling data with models that contain an infinite number of sub-populations 
with different rates of turnover, we have shown that this is in fact feasible. By fitting artificial 
labeling data, we have validated these new models: they generally give good fits to the data and 
converge on average turnover rates that are close to the known average turnover rate. Moreover, 
the new explicit heterogeneity model outperformed the Asymptote model when it came to fitting 
experimental data, especially when the rates of label accumulation and loss are not exponential 
(see Figure [H]). Importantly, due to its relatively general structure, all results of the kinetic 
heterogeneity model are applicable to both non-proliferating and proliferating lymphocytes, 
all having a distribution of turnover values (results not shown). Moreover, because the model 
naturally incorporates the dependence of the rate of label loss on the length of the labeling 
period, this is the first model that can be strictly applied to fit labeling data with different 
labeling periods. 

We have focused our analysis on a particular type of kinetic heterogeneity in which kinetic 
properties of cells of a given subpopulation do not change over time and there is no exchange of 
cells in different subpopulations. However, in some circumstances, such as during an immune 
response, this assumption maybe violated since over the course of infection, lymphocytes do 
changes their kinetic properties over time (e.g., (j2ol)). Taking into account such a type of 
temporal heterogeneity is therefore biologically justified. We have shown that estimates of the 
average turnover rate of the population may depend on the model used (e.g.. Figure S] and [7]B); 
it is therefore important to develop additional (simple) models that take temporal heterogeneity 
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of cell populations into account (12 ll . l22l). We propose that future studies should aim at testing 
multiple models in how well they describe the labeling data and whether these models deliver 
similar estimates of important kinetic parameters such as the average rate of cell turnover. 



4 Methods 



When fitting experimental data, the models were extended to allow for the initial delay in the 
labeling of cells (see also (j^)- For example, including a delay in the Asymptote model (given 
by eqn. ([1])) takes the form 



0, if t < r, 

L{t)={ a(l -e-'^^*-")), if T <t <T + T, (9) 
L(T)e-'^(*-^-^), otherwise. 

To normalize the residuals of the model fits to experimental data, given that the data are 
expressed as proportions, the data and the model predictions were transformed as arcsin(yx) 



where x is the frequency of labeled DNA in the population (1231 ). The models were fitted 
according to the least squares method by using the FindMinimum routine in Mathematica. 
Confidence intervals were calculated by bootstrapping the residuals with 1000 simulations. 
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5 Supplementary Information 

5.1 Derivation of the model with continuous kinetic heterogeneity 

The average label incorporation in a population consisting of n sub-populations is given by eqn. 
([2]). If the number of sub-populations is large {n oo), one could switch from summation to 
integration in eqn. ([2]). If is the fraction of cells in sub-population i with turnover rate dj, 
then, as n — i> oo, = f{d)dd, where the latter is the probability that a randomly chosen cell 
from a population will have a turnover rate in the range {d, d + dd). If Y17=i = 1; then f{d) 
is also normalized to 1. Then for t > T, eqn. ([2]) can be rewritten in a different form 



L{t) = lim 



a,e 



-d^t 



i=l 



POO 

/ fid)e-'^'dd. 
Jo 



(S.l) 



Similarly, from eqn. ([2]), during the delabeling period, the fraction of labeled DNA is given 



by 



POO roo 

L{t) = / /(d)e-'^(*-^)drf - / /(d)e-'^*drf. (S.2) 
Jo Jo 

One can then calculate several important characteristics that determine the change of the 
fraction of labeled DNA over time. First, the initial uplabeling rate can be calculated using 
eqn. (IS. 10 for small t, yielding 



/•oo noo 

L{t) = l- f{d)e-'^'dd^l- f{d){l-dt)dd = dt. (S.3) 
Jo Jo 

where d = df{d)dd, and f{d)dd = 1 by definition. The importance of this result is 
that it demonstrates that for any distribution f{d), the estimated initial rate of uplabeling is 
determined only by the average rate of cell turnover. 

During delabeling, the initial change in the fraction of labeled DNA (for t = T + e with e 
relatively small), can be calculated from: 



L{t) = / f{d)e-''^'-^^dd- / 
Jo Jo 



f{d)e-'^'dd 
15 



f{d){l-de)dd 



/(ci)e-'^^(l - de)dd 

d — df{d)e 



L{T) - eL{T) 



1- / f{d)e 
Jo 



-dT 



dd — e 



d 



df{d)e-''' dd 



l-J,^f{d)e-^^dd\ 



(S.4) 



where L(T) = 1 — f{d)e~^^dd. Then the initial per capita rate of loss of labeled DNA, d*, 
can be calculated for short and long labeling periods. For short labeling periods T — > 0, and 
we find 



^ d- j^df{d)e-''^dd _ d- J,°^dfid)il - dT)dd _d' _ - var(rf) 

1- f{d)e-^^dd l-/o"/(rf)(l-rfr)dd d d ' ^ ' 

where var((i) = d'^ — (d)^ is the variance of the turnover rates in the population. When the 
labeling period is long, T ^ oo, the per capita rate of label loss is 



d-J-df{d)e-^^dd , 

1 - /~ f{d)e-dTdd ' ^ ' 



since all terms e~'^^ ^ as T — oo. This confirms the conjecture of Asquith et al. (jlOl ) that 
d* should approach the average rate of turnover after long labeling period. We now add that 
the maximal difference between the average turnover rate, d, and the loss rate of labeled cells, 
d*, is set by the distribution of turnover rates in the population. 



5.2 Particular solutions of the kinetic heterogeneity model 

For several simple distributions of the turnover rates, we can obtain analytical solutions for the 
change in the fraction of labeled DNA during the labeling experiment. 

Exponential distribution. In this case f{d) = {l/d)e-'^/'^, where d is the average turnover 
rate in the population. Using eqn. (IS.ip and (lS.2p . we find 



L(i) 



1-1 

d 



-dt 



-"^/^dd = 1 



1 



{t + l/d)d 



l + dt' - ' 



(S.7) 
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1 

d Jo 

L{T) 



d Jo 



■'^/'^dd 



1 



1 



l + d{t-T) 1 + dt 



1 + dT 



[l + dt){l + d{t-T)) 



t > T. 



(S. 



When dt -C 1, the fraction of labeled DNA is simply L{t) ~ dt, i.e., the initial rate of increase 
is again given by the average turnover rate d. The initial rate of decline during delabeling is 
less clear. Let us define t = T + e where pe <^ 1. Then after cessation of label administration, 
using Taylor's expansion we find 



Lit) 



1 



1 



1 — de 



l + de l + dT + de 1 + dT 

L{T) - L{T)d (l + e + o{e) 



+ 



de 



[l + dTf 



+ o(e) 



(S.9) 



where L{T) = 1 — 1/(1 + dT). This expression shows that the initial per capita rate of loss of 



labeled DNA, d* = dyl + r^j' decreases with increasing length of the labeling period (llOl ). 
If dT <^ 1 (short labeling), the initial per capita decay rate is d* ~ 2d. 

Gamma distribution. In this case f{d) = X{Xd)'^~^e~'^'^/{k — 1)!, where A and k are 
the scale and shape parameters, respectively, and the average rate of cell turnover d = k/X, 
cr^ = d'^/k, and CV = ad/ d = 1/ \fk. Note that when = 1, the gamma distribution is identical 
to the exponential distribution. Substituting the gamma distribution in eqn. (IS. Ill , we find 



1 - 
1 - 
1 - 



A* 



(^-l)!7o 
A'^ 



{k-m + XY }o 
X^ 



d^-\t + Xfe-'^'^'+^^dd 



it + A)^ 



1 - 



1 + 



dt 
1 



n -k 



t < T. 



(S.IO) 



Using eqn. flS.2p and proceeding similarly, we find the change in the fraction of labeled DNA 
during delabeling: 



L{t) 



1 + 



d{t - T) 
k 



1 


' dt' 







t > T. 



(S.ll) 
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(S.12) 



The initial rate of increase in the fraction of labeled cells is also independent of the length 
of the labeling period and, initially, for t — e (such as de <^ 1) is 



Lie) = 1 - 



1 + 



de 



1 - (1 - kde/k) + o{e) = de + o{e). 



(S.13) 



as expected. The initial per capita rate of loss of labeled DNA is somewhat more complex. For 
times t — T + e such as de <^1, using Taylor's expansion, we find 



de 



1 + 



d{T + e) 
k 



-k 



1 - 



k) 



— de 



1 - 



(S.14) 



It is useful to rewrite this expression in terms of L{T): 



L{t) = L{T)-L{T): 



de 



dT/k 



[1+dT/k) 



k+l 



(1 + dT/k)'' - 1 



L{T){l-d*e). (S.15) 



where d* 



l+dT/k 



k + l_ 



(1+dT/k) 
(l+dT/fc)*-! 



is the initial per capita loss of labeled DNA. For short labeling 

{dT ^ 1), the initial per capita decay rate is d* ^ d{k + l)/k, and as the shape parameter 
k becomes larger, the decline rate d* approaches the average proliferation rate d. For long 
labeling periods (dT 1), d* ^ d, as expected. 

Truncated gamma distribution. Under some circumstances the distribution of turnover 

rates may allow for too high rates of cell turnover. To circumvent this problem one may use 
a truncated distribution. For a gamma distribution truncated at maximal value dmax, the distri- 
bution is similar as above with an added normalization constant C, f{d) = C~^A(Aci)'^^^e^^°'/ (A;— 
1)!. The constant C is found by normalizing the probability distribution 



C= / f{d)M 
Jo 



1 - 



r(/c, d^nax^) 

{k-iy. ' 



(S.16) 
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where T{k, d) = ^dx is an incomplete gamma function. The average turnover rate 

then has to be calculated numerically 



d = C'^ / df{d)dd 
Jo 



X 



k\ 
(k-l)\ 



For the fraction of labeled nucleotides, we proceed as in eqn. (IS.lOp and obtain 



(S.17) 



L{t) 



'k-l-dt-Xd 



d'^-'e 



C{k-l)\J, 

yk rdmax{t+X) 

c{k-i)\{t + xY j, 
x'^ 



dd 



x^ ^e "^dx 



C{k-l)\{t + XY 



{k-l)\- / x^^-V^dx 



X' 



{t + xy 



X 



ik-l)\-Tik,dmaA^ + t)) 
{k-l)\-T{k,dma.X) 



t < T. 



During delabeling {t > T), we proceed as in eqn. fIS.lip and find 



A* 



{x + t-ry 



{k-l)\-T{k,dmaAX + t-T)) 
{k - 1)! - T{k,dmax>^) 



x'^ 



it + xy 



X 



ik-iy-T{k,d^a,{X + t)) 
{k - ly. -T{k,dmax>^) 



(S.19) 



Since lim^^^^^^oo r(A;, 



, at dn 



CXD, 



the fraction of labeled nucleotides becomes identical to eqn. (jlj) with A = k/d. 

5.3 Non- parametric estimates of the distribution of proliferation 
rates in the population 



Mathematically, the last term in eqn. (IS.ip is a Laplace transformation of f{p). This result 
stems from the assumption of exponentially distributed inter-division times of cells and raises 
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the intriguing possibility that from the change in the fraction of labeled DNA during label 
administration, one can estimate the distribution of proliferation rates in the population. Let 
us denote as the Laplace transformation of f{d): 



nOO 

nt) = C[fid)] = / f{d)e-'''dd. 
Jo 



(S.20) 



From eqn. ( IS.ll) . one finds the distribution of proliferation rates in the population using the 
inverse Laplace transformation: 



f{d) = C-'[l-Lit)]. (S.21) 

where L{t) is a curve describing the change of the fraction of labeled DNA during label ad- 
ministration. Such a curve could be obtained in several ways, for example, by interpolating 
the data. Application of this method for the analysis of experimental data will be published 
elsewhere. 
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Gamma model 



Data generated using: 

Exponential model 



Two-populations model 





7 days 


15 days 


7 days 


15 days 


7 days 


15 days 


d 


0.076 


0.068 


0.082 


0.08 


0.04 


0.039 




0.068—0.085 


0.06—0.077 


0.073—0.096 


0.07—0.092 


0.034—0.046 


0.032—0.047 


d 


0.121 


0.117 


0.203 


0.199 


0.052 


0.054 




0.105—0.14 


0.099—0.137 


0.17—0.238 


0.169—0.231 


0.04—0.066 


0.038—0.071 


d 


0.09 


0.09 


0.11 


0.12 


0.05 


0.05 




0.082—0.104 


0.076—0.1 


0.105—0.122 


0.104—0.137 


0.042—0.055 


0.040—0.055 


a 


0.83 


0.76 


0.51 


0.48 


1 


1 




0.773—0.908 


0.70—0.82 


0.488—0.536 


0.456—0.51 


0.857—1.0 


0.835-1.0 


d 


0.097 


0.099 


0.154 


0.197 


0.044 


0.044 




0.085—0.112 


0.085—0.117 


0.14—0.169 


0.152—0.271 


0.038—0.052 


0.036—0.054 


k 


0.594 


0.44 


0.197 


0.159 


1.329 


1.132 




0.478—0.768 


0.367—0.536 


0.184—0.211 


0.134—0.183 


0.796—2.947 


0.679—2.901 



J3 

4J 



T3 

EC 

Q 



CD 

o 
ft 

a 

IB 



CD 

o 
ft 

X 



o 



Table 1: Estimates of the parameters after fitting three models to three sets of artificial data. Rows 
correspond to different models used to fit the data, and columns correspond to the models used to 
generate the data. We show parameter estimates of the Asymptote model, the Exponential model 
(in which a fraction a of cells have exponentially distributed turnover rates), and the Gamma model 
(with gamma distributed turnover rates). Data were generated using the Gamma model (eqn. dH)), 
the Exponential model (eqn. ([8])), and the Two-populations model (eqn. ([2])). The 95% confidence 
intervals, that are shown below the mean values, were obtained by bootstrapping the residuals with 
1000 simulations. Data have been generated for two labeling periods, of 7 and 15 days, respectively. 
For all data, the average rate of turnover was fixed at 0.1/day. Other parameters used to generate the 
data are: k = 0.5 (Gamma model), a = 0.5 and da = 0.2/day (Exponential model), di = 1/day, and 
a = 0.07 (Two-populations model). 
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Data fitted with: 




AsviTiDtotp rnndpl 


Exponential model 


Ga-mnici model 


dldr.lk 


2.87 (2.22 — 3.65) 


1.94 (1.46 — 2.6) 


0.12 (0.09 — 0.17) 


di, % day^^ 


0.57 (0.5—0.67) 


0.59 (0.51—0.7) 


0.62(0.53—0.78) 


d2 


0.41 (0.34—0.5) 


0.44 (0.36—0.53) 


0.43 (0.34—0.54) 


ds 


0.38 (0.31—0.47) 


0.41 (0.33—0.5) 


0.37 (0.27—0.48) 


di 


0.41 (0.33—0.5) 


0.44 (0.36—0.54) 


0.43 (0.32—0.57) 


n, day 


1. (0.91—1.51) 


1. (0.93—1.53) 


1. (0.93—1.57) 




0.78 (0.28—1.) 


0.81 (0.35—1.) 


0.8 (0.31—1.) 




1.97 (1.— 2.65) 


2.06 (1.14—2.7) 


1.87 (0.73—2.67) 




1.7 (0.98—2.45) 


1.83 (1.-2.54) 


1.79 (1.-2.58) 


RSS, 10^3 


6.19 


5.94 


5.87 



Table 2: Average turnover rates of 004"*" T cells from four healthy humans as estimated by fitting 
the data from Mohri et al. (0) using the Asymptote model, the Exponential model, and the Gamma 
model. The best fits of the models resulted in different average rates of cell turnover di and initial delays 
of labeling r^. Other parameters, that could be assumed to be identical between different patients, 
are the death rate of labeled cells d (Asymptote model), the rate of turnover da of the turning-over 
sub-population in the Exponential model, and the shape parameter k in the model with gamma 
distributed turnover rates. For the model with gamma distributed turnover rates, an asymptote level 
a = 1 provided the best fit of the data. The quality of the fit is illustrated by the residual sum of 
squares (RSS). The 95% confidence intervals were obtained by bootstrapping the residuals with 1000 
simulations. 
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Data fitted with: 

Asymptote model Exponential model Gamma model 


ai/a/ki 
a2/a/k2 


0.08 (0.07—0.10) 
0.13 (0.09—0.73) 
0.26 (0.10—1.0) 
0.10 (0.07—0.28) 


0.13 (0.12—0.16) 
0.13 (0.12—0.16) 
0.13 (0.12—0.16) 
0.13 (0.12—0.16) 


0.033 (0.028—0.039) 
0.116 (0.066—0.341) 
0.347 (0.099—10^) 
0.082 (0.05—0.155) 


di, % per day 

d2 

d3 

di 


0.36 (0.31—0.43) 
0.23 (0.18—0.30) 
0.21 (0.17—0.31) 
0.22 (0.18—0.3) 


0.48 (0.39—0.57) 
0.27 (0.21—0.34) 
0.29 (0.22—0.39) 
0.24 (0.2—0.32) 


0.62 (0.53—0.71) 
0.23 (0.19—0.28) 
0.2 (0.17—0.27) 
0.23 (0.19—0.28) 


Tl 
T2 
T3 

ta 


0.99 (0.73—1.38) 
0.65 (0.-0.94) 
0.85 (0.— 1.99) 
0. (0.— 0.63) 


1.76 (1.33—1.95) 
0.76 (0.11—0.99) 
1.73 (0.59—2.4) 
0. (0.— 0.64) 


1.89 (1.75—1.98) 
0.65 (0.28—0.92) 
0.82 (0.-1.81) 
0. (0.-0.55) 


RSS, 10-3 


3.4 


3.85 


1.56 



Table 3: Average turnover rates of 008"*" T cells from four healthy humans as estimated by fitting 
the data from Mohri et al. ([3) using the Asymptote model, the Exponential model and the Gamma 
model. The best fits of the models resulted in different parameter estimates for all patients, with the 
exception of the fraction of turning over cells a in the Exponential model (which was fitted as one 
parameter for all patients). As for 004"*" T cells, in the model with gamma distributed turnover rates, 
the asymptote level a = 1 provided the best fit of the data. The shown 95% confidence intervals were 
obtained by bootstrapping the residuals with 1000 simulations. 
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Figure 1: A cartoon of the model with explicit kinetic heterogeneity. In the model, the population 
of cells consists of n sub-populations with different rates of turnover. In the i*'* sub-population, there 
is a source of new cells that enter the cell population at rate Si cells per day, cells divide at rate pi per 
day, and die at rate di per day. To maintain the size of all sub-populations constant, di — pi = si/Xi 
for every sub-population i, where Xi is the number of cells in the i*^ sub-population. In this model 
we assume that the source produces only labeled cells during the labeling phase, and delabeled cells 
during the unlabeling phase (jlll ). 
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Figure 2: Model predictions for exponentially distributed turnover rates. We have plotted the 
changes in the fraction of labeled DNA according to the explicit kinetic heterogeneity model with 
exponentially distributed turnover rates (eqn. ([5]), mean d = 0.1/day). Predicted changes are shown 
for a short labeling period (T = 1 day, solid line) and a long labeling period (T = 30 days, dashed line) 
on a linear (panel A) and a logarithmic (panel B) scale. The initial uplabeling rate is independent of 
the length of the labeling period and is given by d. The initial rate of delabeling, in contrast, depends 
on the length of the labeling period and is approximately twice as fast in the case of short-term labeling 
as compared to long-term labeling. 



24 




Figure 3: Model predictions for gamma-distributed turnover rates. We have plotted the changes 
in the fraction of labeled DNA according to the kinetic heterogeneity model with gamma-distributed 
turnover rates (eqn. ^) with average turnover rate d = 0.1/day on a linear (panel A) or logarithmic 
(panel B) scale. Predicted changes are shown for different values of the shape parameter k. Larger 
values of k correspond to a more symmetric distribution (Panel C). For low values of the shape 
parameter k, the loss of labeled DNA after label cessation is biphasic, which is most clearly visible 
on a logarithmic scale for A; < 1 (panel B). This characteristic of the kinetic heterogeneity model 
differs from the Asymptote models which have a constant per capita rate at which labeled DNA is 
lost. Note that for shape parameters k < 1, the distribution of turnover rates f{d) becomes extremely 
skewed with most cells undergoing hardly any division and relatively few cells undergoing extremely 
many rounds of division (panel C). Panel D gives the cumulative contribution of sub-populations with 
a particular turnover rate d to the average rate of turnover of the population d. The vertical line 
shows the value of the average proliferation rate d. For high values of the shape parameter (k = 5), 
the cell sub-populations with turnover rates that are somewhat lower or higher than d give the main 
contribution to the average turnover rate. In contrast, for low values of k {k = 0.01), the major 
contribution to the average turnover rate comes from sub-populations with extremely rapid turnover 
rates {d^ d); about 50% of the average turnover is due to a few sub-populations with turnover rates 
that exceed 10 per day, which is biologically unrealistic. 
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Figure 4: Fitting artificial data (black dots) with the Asymptote model (eqn. ([T|), solid black lines), 
the Exponential model, in which a fraction a of the cells have exponentially distributed turnover rates 
(eqn. ([8]), small red dashed lines), and the Gamma model with gamma distributed turnover rates 
(eqn. ([4]), large green dashed lines). Data were generated using the Gamma model (panel A&B), the 
Exponential model (panel C&D) and the Two-populations model (panel E&F, eqn. (H])), respectively. 
Thin blue lines show the exact curves of the models that were used to generate the data. The different 
models were fitted to 11 datapoints taken from these predicted curves after having added noise to these 
data points. Noise was added by a relative change of the predicted value with a normally distributed 
error (with standard deviation of the distribution a = 0.1). The models were fitted to data from 
artificial labeling experiments in which the label was administered for 7 (left panels) or 15 (right 
panels) days. Parameter estimates providing the best fit are shown in Table [H and the corresponding 
estimates of the average rates of cell turnover d are shown in Figure O Parameters used to generate 
the data are also given in Table [H ofi 
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Figure 5: By fitting the artificial data described in the text, we estimated the average turnover rate 
using three models: the Asymptote model (eqn. ([I|)), the Exponential model (in which a fraction 
a of cells have exponentially distributed turnover rates, eqn. ([8])), and the Gamma model (with 
gamma-distributed turnover rates, eqn. (131) )• Estimated mean values and 95% confidence intervals 
obtained by bootstrapping the residuals with 1000 simulations are shown. Data were generated using 
the Gamma model (□), the Exponential model ( ) and the Two-populations model (■). Labeling 
periods were 7 (panel A) and 15 (panel B) days. Horizontal dashed lines denote the actual average 
rate of lymphocyte turnover in all data, d = 0.1/day. Note that in this example, the Asymptote 
model always underestimated the average rate of cell turnover, and that there is a systematic 2-fold 
underestimation of the average turnover by all models when the data from the Two-populations model 
were fitted. This is because all three models fail to describe the relatively rapid accumulation of the 
label at early time points (see Figure [3E~F) . 
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Figure 6: Fitting the deuterium labeling data of CD4^ (top rows) and CD8"^ (bottom rows) T cells 
in four healthy humans with three models: the Asymptote model (panels A and D), the Exponential 
model, in which a fraction a of cells have exponentially-distributed turnover rates (eqn. (I8|), panels B 
and E), and the Gamma model, with gamma-distributed turnover rates (eqn. (jlj), panels C and F). 
Experimental data obtained from Mohri et al. (0) are shown as symbols and the curves are the best 
model fits. The sum of squared residuals of the model fits to the data on the dynamics of CD4'*' T 
cells are (6.19,5.94,5.87) x 10~^ for the Asymptote model, the Exponential model and the Gamma 
model, respectively. The sum of squared residuals of the model fits to the data on the dynamics 
of CD8^ T cells are (3.4,3.85,1.56) x 10^'^ for the Asymptote model, the Exponential model and 
the Gamma model, respectively. Note that the two explicit kinetic heterogeneity models describe 
these data with similar (Exponential model) or even better (Gamma model) quality compared to the 
Asymptote model. 
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Figure 7: Estimates of the average turnover rates of 004"*" (panel A) and CDS'*" (panel B) T cells 
in four healthy humans obtained by fitting with three different models: the Asymptote model (□, the 
Exponential model (-,), in which a fraction a of the cells have exponentially-distributed turnover rates, 
and the Gamma model (■) with gamma-distributed turnover rates. Best fits of the data are shown 
in Figure [6l and estimates of all parameters of the models are shown in Appendix (Table [2] and [3]). 
Confidence intervals were obtained by bootstrapping the residuals with 1000 simulations. Horizontal 
dashed lines denote the mean of the estimated turnover rates in all patients and all models for CD4"'" 
{d = 0.46% per day) and CD8^ {d = 0.29% per day) T cells. Note that all models deliver very similar 
estimates for the average turnover rate d, with the exception of the estimated CD8+ T-cell turnover 
rates in individual cl which are highly mo del- dependent. 
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